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A population evolving in an inhomogeneous environment will adapt differently to different regions. 
We study the conditions under which such a population can maintain adaptations to a particular 
region when that region is not stationary, but can move. In particular, we study a quasispecies 
living near a favorable patch ("oasis") in the middle of a large "desert." The population has two 
genetic states, one of which which conveys a relative advantage while in the oasis at the cost of 
a disadvantage in the desert. We consider the population dynamics when the oasis is moving, or 
equivalently some form of "wind" is blowing the population away from the oasis. We find that 
the ratio of the two types of individuals exhibits sharp transitions at particular oasis velocities. 
We calculate an extinction velocity, and a switching velocity above which the dominance switches 
from the oasis-adapted genotype to the desert-adapted one. This switching velocity is analagous to 
the quasispecies mutational error threshold. Above this velocity, the population cannot maintain 
adaptations to the properties of the oasis. 



Spatial inhomogeneities in the environment are often 
essential to understanding the dynamics of natural pop- 
ulations. Populations may for example be confined to 
limited reserves or live in an environment with gradients 
in resources or large-scale inhomogeneities in habitability. 
The evolution of a population in an inhomogeneous en- 
vironment is particularly interesting. When individuals 
move over the environment sufficiently rapidly relative to 
the length scale of the inhomogeneities, they may all see 
and adapt to an "averaged" environment. When this is 
not true, however, the dynamics can be much more com- 
plex. Some individuals may randomly see primarily one 
part of the range while others see other parts. The pop- 
ulation in the distant future will likely be dominated by 
the descendents of a few exceptionally lucky individuals 
in the present, who will have seen an unusually favorable 
subset of the set of possible environments. Thus it is not 
at all clear a priori exactly how different regions of the 
environment will influence the evolution. 

In this paper, we consider a simple environment with 
two regions, a favorable "oasis" in a large less favorable 
(or unfavorable) "desert." The favorable area could be 
realized as a game reserve, a region of favorable climac- 
tic conditions, or a patch of light, among other things. 
Provided that the favorable region is stationary and suffi- 
ciently large, a population will adapt to the special prop- 
erties of this region. However, in many cases the oasis 
will move at some typical speed v, for example due to sea- 
sonal weather patterns or human intervention. (Equiva- 
lently, the population could be blown across the oasis at 
speed v by some form of "wind"). In this situation, there 
is a velocity beyond which the population will not be 
able to maintain adaptations to the properties of the oa- 
sis. Rather, the evolution will be dominated only by the 
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desert environment. We focus in this paper on calculat- 
ing this velocity, which we call the "switching velocity." 

Recent work has considered the dynamics, though not 
the evolution, of a population being blown by some form 
of convective "wind" across a nonuniform environment. 
In these models, the equation for the population density 
at a point (x, t) is 

=DW 2 c-v-Wc + r(x)c, (1) 

at 

where D is the diffusion constant, v is the drift veloc- 
ity, and r(x) is the growth rate. This equation thus 
describes a population multiplying in response to some 
spatially heterogeneous environment while diffusing and 
drifting. It can be analyzed using techniques adapted 
from an analysis of non-Hcrmitian Schrocdinger-like op- 
erators 0,0 ■ A nonlinear saturation term can be added, 
and is important for some purposes. 

Ref. 0] examines this model in two dimensions in the 
case where r(x) is random, with only short-range cor- 
relations. In the limit of small v (v -C vf = 2>/aD, 
where a is the average of r(x)), the population is domi- 
nated by a few colonies that grow up around "hot spots," 
regions which happen to have higher growth rates than 
surrounding areas. As v increases, individual colonies on 
less prosperous hot spots are blown away in a series of 
"delocalization" transitions. In the limit of large v all 
organisms are blown away from individual hot spots and 
convect across the environment. At long times, the pop- 
ulation is dominated by those individuals that happened 
to take very special paths through the random environ- 
ment, travelling through a disproportionate number of 
the hot spots. Thus this is an example of a system in 
which a typical individual will evolve in response to a 
very special subset of the environment. 

In subsequent work, Dahmcn et. al. examined the 
transition between the small and large v regimes by look- 
ing at a model of a sing le "hot spot" 0. These authors 
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took the growth rate r(x) to be large within a small hot 
spot, or "oasis," and smaller (and possibly negative) in 
the surrounding "desert." They analyzed the dynamics 
near the extinction and derealization velocities, where 
the population is blown off of the oasis. Their predictions 
have been qualitatively confirmed by recent experiments 
on Bacillus subtilis growth Q. 

In this paper, we extend this work to consider the dy- 
namics of a population of multiple types of individuals in 
a desert-oasis environment. In particular, we consider a 
quasispecies model with two types of individuals. Thus 
c(x, t) becomes a two- vector 



c(x, t) 



ci{x,t) 
c 2 (£, t) 



(2) 



and r{x) becomes a two-by-two matrix with diagonal el- 
ements specifying the growth rates of the two genotypes 
and off-diagonal elements specifying the mutation and 
back-mutation rates. Genotype 1 represents the popula- 
tion of some "ideal" genome sequence, while type 2 rep- 
resents all other less ideal sequences. There is mutation 
back and forth between the two, typically with a muta- 
tion rate away from the ideal sequence greater than the 
mutation rate towards it. Without any spatial variation, 
(i.e. with r(x) independent of x), this is a simple quasis- 
pecies model. A few recent papers have examined other 
versions of a spatial quasispecies model, with many differ- 
ent possible types of individuals representing all possible 
Hamming distances from the ideal sequence |6|, [jj . Here, 
we focus on a simple two-type spatial model, although it 
is straightforward to generalize our results. 

We assume that the population c\ has some gene or se- 
quence which conveys an advantage in the environment 
of the oasis, but has an overall cost which makes it dis- 
advantageous in the desert. The population c 2 consists 
of all individuals who have lost the function of the gene 
by one or more mutations. Individuals of type 1 are thus 
those that have adapted to the special properties of the 
oasis, while type 2 individuals are desert-adapted. 

We examine the population dynamics as we change the 
velocity v. We find two important transitions. When the 
growth rate in the desert is negative, there is an extinc- 
tion velocity v e where the entire population goes extinct. 
When the genotype 2 desert growth rate is positive, there 
is a "switching" velocity v s where the behavior of the 
population changes dramatically. Below this velocity the 
population (particularly within the oasis) is dominated 
by genotype 1. Above it, the oasis- adapted genotype is 
outcompeted, genotype 2 dominates, and C1/C2 — > at 
long times. This transition is a velocity-driven analog 
of the classical mutational error threshold in non-spatial 
quasispecies models 0, 0, ^| . 

These results have interesting implications. If the oasis 
represents some part of a species' habitat, the switching 
velocity is simply how quickly this can move before the 
species can no longer maintain adaptations to the prop- 



erties of this aspect of its range. Our analysis is also a 
first step towards understanding the spatial quasispecies 
model in a random environment, where the population 
does not evolve in response to the averaged environment 
but rather in response to a different, highly selective sub- 
set of the environment. 

The outline of this paper is as follows. In section U we 
give a detailed description of our model. In section ITT1 we 
outline the calculation of the critical velocities and the 
behavior of the populations in the different regimes. In 
section UTTl we compare our analytical results with com- 
puter simulations. Finally, in section llVl we discuss the 
main biological implications of these calculations. 



I. MODEL 

We consider a population with two types of individuals 
whose densities are described by 



c(x, t) 



ci(x,t) 
c 2 (x,t) 



(3) 



The dynamics is exponential growth (or decay) with dif- 
fusion and convection 



dc(x, t) 
dt 



DW 2 c-v-Wc + TZ(x)c, 



(4) 



where 72.(x) is a two-by-two matrix. We are primarily 
concerned with the extinction and derealization transi- 
tions. Thus for the bulk of our analysis, we neglect the 
possibility of a nonlinear saturation term. In scction liV Al 
we discuss the consequences of such a nonlinearity. 

We define the growth and death rates of individuals 
of type i within the oasis to be on and respectively. 
Outside the oasis, the growth and death rates are defined 
to be f3i and Si. The mutation rate from type 1 to type 
2 is defined to be /ii, and the back-mutation rate is // 2 . 
Based on these definitions, the matrix IZ(x) is given by 



7l{x) 



osi — 71 - A*i"i 



A*2«2 
Ct2 — 72 — 



(5) 



inside the oasis and similarly (with a — > (3 and 7 — > 5) in 
the desert. 

For simplicity, we assume that the growth rates within 
the desert and the oasis are the same, i.e. an = Pi. The 
desert is less hospitable because of a higher death rate, 
not a lower growth rate. This simplifies the analysis but 
is not a serious limitation; the results for a 7^ (3 are 
straightfoward to calculate by the same methods. We 
next define = a 1 — 71 and bi = Pi — Si. For convenience 
we also assume that a± — a% and et 2 = a 2 (i.e. 71 = 7 2 = 
0). It is easy to relax this assumption, but the results 
more transparent with it in place. These simplifications 
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yield the matrix 



h - Mi«i M202 



(inside oasis) 



in the oasis and similarly in the desert, 



(6) 



n{x) = 



h - l-Ho-i 



b 2 - fl2d2 



(outside oasis). 



(7) 

We assume that a± > &i and a 2 > b 2 , so that we have a 
beneficial oasis and a harmful desert. We further assume 
that a\ > a 2 and b 2 > b\ (which also implies a 2 > b\) so 
that genotype 1 is better in the oasis and genotype 2 is 
better in the desert. Absent these inequalities, one or the 
other type will unequivocally dominate the population 
(up to the mutational error threshold) independent of the 
drift velocity. In that case, we can ignore the inferior type 
at long times and the problem reduces to that studied 
in Q. Our assumptions imply that individuals of type 
1 have some function that conveys an advantage inside 
the oasis, at the cost of a disadvantage elsewhere. Our 
analysis will determine whether or not this function can 
be maintained in the face of mutational pressure when 
the oasis is moving at velocity v. 

We could of course make any number of assumptions 
about the geometry of the desert and the oasis. We will 
consider for simplicity a one-dimensional system with an 
oasis of width W in the middle of an infinite desert, as 
depicted in figure la. This analysis also describes the 
long time behavior of the two-dimensional system with 
geometry as described by figure lb. If the initial condi- 
tions are uniform in the y direction, the one-dimensional 
system of figure la is equivalent to the two-dimensional 
system of figure lb. If the conditions are nonuniform, 
the two become equivalent after a time of order £ 2 /D £|. 
More complicated geometries are certainly possible, and 
we discuss these briefly in section ITV BI 

We assume throughout that differential equations are 
an adequate representation of these biological processes. 
We neglect discreteness in population number, which 
may have some importance near the extinction transi- 
tion. The effects of discreteness may be analyzed using 
the methods of 



s oi discreteness may DC ; 

in ei in in e ana 



II. CALCULATION OF THE POPULATION 
DYNAMICS 

We analyze our model using the non-Hermitian meth- 
ods of Refs. 0,H 13- We first rewrite our system in the 
form 



dc(x, t) 

at 



Cc, 



(8) 



where C is given by £ = C ln 9(W — \x\) + C QUt 6(\x\ — W), 
and W is the width of the oasis. The function 6{y) = 1 
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FIG. 1: (a) shows the one-dimensional system we consider, 
for the case of a deadly desert. The dashed line is the growth 
rate of genotype 1, the dotted line the growth rate of genotype 
2. This is an approximation to the two-dimensional system 
shown in (b), valid when the initial conditions are uniform in 
the y direction or the time is greater than £ 2 /D. 



if y > and otherwise, and 

Dd 2 — vd x + ai — /iiai 



/ii<xi Dd x 

DO 2 — vd x + bi — (iidi 

/iiai Dd^ 



- vd x + a>2 — fl2d2 

vdx + 62 — M2C12 



C is non-Hermitian, but we can diagonalize it with a 
system of left and right eigenfunctions 4>n(x) and <p^(x), 
with their (common) eigenvalues T n . These eigenfunc- 
tions (which we also call "states") satisfy the orthonor- 
mality condition 



0m 0*0 • <l>n( x ) dx = S " 



(9) 



Using this result, we can write the initial condition as 
a linear superposition of the right eigenfunctions, i.e. 
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c(x, t = Q) = J2 n c n<t>n( x )i where 



Thus for v > 0, 



c n = I d d xcj)^(x) ■ c(x, t = 0). 



(10) 



a vx/2D-K n \x\ 



(15) 



We can then immediately write down the solution valid 
for all times, namely 



c(x,t) = Vc„</>^(x)e 



(11) 



This result has a clear biological interpretation. Ini- 
tially, the population is in some particular arrangement, 
expressible as a linear combination of the different eigen- 
functions. This combination will necessarily include the 
eigenfunction corresponding to the largest T n = T gs 
(which is always real) [l^. This eigenfunction is special, 
and we refer to it as the "ground state" As time 

passes, the population distribution looks more and more 
like the ground state. This distribution will grow (or die 
out) exponentially with rate r gs . The other </>^j may be 
important initially, but since they grow more slowly than 
the ground state, they soon become irrelevant. Thus un- 
derstanding the ground state function and its eigenvalue 
are the key to understanding the long time behavior of 
the population. 

It remains to solve for the eigenfunctions and eigen- 



When v < 2Dn n , the eigenfunctions vanish at infinity, 
as they should. However, for v > 2Dn ni this func- 
tion blows up at infinity, which is unreasonable. The 
correct eigenfunctions have a different character when 
v > 2Dn n = u*. 

Ref. pj shows that the transition at v* is a delocal- 
ization of the corresponding eigenfunction. For v < i;*, 
the eigenfunctions are localized around the oasis, but for 
v > i>* they are delocalized. This makes intuitive sense. 
For small velocities, the population will tend to cluster 
around the oasis, but for larger wind velocities it gets 
blown off the oasis and must live by drifting across the 
desert. The behavior of the eigenvalues r n near this delo- 
calization transition is striking. Up to u*, the eigenvalue 
is simply equal to T^ _u — j^, but beyond this point 
this relation no longer holds. Instead, T n jumps off the 
real axis at u* , becoming complex, and the eigenfunctions 
become broad delocalized states extending through the 
desert 0, 0, Q . We denote the value of T n at which this 
occurs by T* . As we continue to increase v above u* , the 
real part of r n stays approximately constant, although 
the imaginary part does change. From the gauge trans- 



values 0£(ar),0*(a:), and r n . For the case v = 0, the formation relationship, we have T* = F 
solution is straightforward and will be discussed in de- 
tail below. For nonzero v, we make use of the fact that 
the eigenfunctions of C v= q are related to the eigenfunc- 
tions of C by an "imaginary gauge transformation" 
That is, if 4>n v =o( x ) ^ s a right eigenfunction of C v= o with 
eigenvalue T„, then 

4>* v (x) = e°*' 2D <f>% v=0 (x) (12) 

is a right eigenfunction of C, with eigenvalue 

K = K =0 -^, (13) 

as can be verified by allowing C to act on <j>^. Similar 
expressions hold for </>^ v . The eigenvalues r^ =0 are all 
real. Thus, eigenfunctions for v > are very similar 
to the eigenfunctions for v = Q. The genotype 1 versus 
genotype 2 composition of the states is not altered at 
all. The only change is that the wind causes a distortion 
of the population in the direction of the wind, and the 
growth rates of the states shift downward "rigidly" (i.e. 
independent of n) by an amount ^ . 

However, this procedure works only for small v. To 
see this, consider the behavior of the eigenfunctions as 
x — > oo. We expect (and will soon verify) that the v — 
eigenfunctions far from the oasis decay exponentially, 



4D 



,-K n \x\ 



(14) 



How- 
ever, as we will see, the structure of the n-dependence of 
v n and r„ is such that there are only two different values 
of r*. This is a crucial point. In our problem, the states 
will divide into those dominated by genotype 1 and those 
dominated by genotype 2. As we will show, states domi- 
nated by genotype 1 delocalize at T\ = (n), the spatial 
average growth rate of the first genotype, which up to 
finite size effects is just b\. = (r 2 ) ~ b 2 plays the 
same role for states dominated by the second genotype. 
By our assumptions about the parameters, TZ, > r*. An 
example of eigenvalue spectra for several values of v is 
given in figure [21 

Each localized state has some particular r^ =0 , and for 
nonzero v its eigenvalue becomes T n = T^ =0 — fjj- As 
v increases, T decreases until the state delocalizes, at 
r\ for those states dominated by type 1 and at T 2 for 
those states dominated by type 2. We will see that the 
ground state for v = is dominated by genotype 1. As 
v increases there comes a critical point when this ground 
state eigenvalue T gs crosses T^. Beyond this point it 
is no longer the ground state. Rather, the delocalized 
genotype-2 dominated state with eigenvalue T 2 has the 
highest r, and this state determines the population dy- 
namics. We call this critical velocity the "switching ve- 
locity," where the dominance switches from genotype 1 
to genotype 2. This switch is a type of quasispecies tran- 
sition, caused not by exceeding a mutation rate error 
threshold, but by exceeding a critical velocity. This rea- 
soning is illustrated in figure [3] 
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Typical Eigenvalue Spectra 
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FIG. 2: An example of eigenvalue spectra, shown here for 
3 values of v = — = — t= with a\ = l,b\ = — l,a,2 

0.6,6 2 = -0.6, ^1 = 0.01,^2 = 0.001, and D = 0.25, where 
the overbars indicate the non-dimensionalized parameters dis- 
cussed in Appendix C. To allow easy distinction from the 
v = results, the localized v — 0.9 eigenvalues have been 
shifted slightly upwards. Note the rigid shift of the localized 
eigenvalues (those on the real axis) to the left as we increase 
v. This is highlighted by the three v — states marked by 
upwards arrows, which for v = 0.9 are shifted into the three 
states marked by downards arrows. The two (and only two) 
derealization transitions T* ~ — 1 and T2 ~ —0.6 are also 
clearly visible. Note that these transition points, and all the 
delocalized states, do not shift to the left as v increases. On 
the left half of the spectrum (omitted from the figure), the 
delocalized eigenstates form closed loops, an artifact of the 
computational discretization used to produce this figure. 



If the growth rate for genotype 2 in the desert is neg- 
ative, then Ft, < 0. The switching behavior will then be 
difficult to observe experimentally, as both genotypes will 
be going extinct when it happens. In this case, the bio- 
logically interesting transition occurs at the "extinction" 
velocity where the growth rate of the ground state passes 
through 0. Beyond this velocity, all the states have nega- 
tive growth rate, so neither genotype can survive. When 
the growth rate for genotype 2 in the desert (& 2 ) is posi- 
tive, the switching velocity becomes biologically relevant. 
In this case, there is no extinction velocity. This is be- 
cause the delocalized states do not shift to lower T as 
v increases, and at least one genotype-2 dominated de- 
localized state has eigenvalue T% ss 6 2 > 0. Intuitively, 
this makes sense because the population, dominated by 
genotype 2, can survive in the desert at arbitrarily large 
velocities. 



Typical Eigenvalue Spectrum 
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FIG. 3: A single eigenvalue spectrum, with the genotype-1 
and genotype-2 dominated states distinguished. This is just 
the v = 0.9 spectrum from figure above. Note the state 
with the largest Re(r) (i.e. the ground state) is a genotype-1 
dominated state. This v is just below the extinction velocity, 
so T gs is just above the extinction threshhold, indicated by 
the dashed vertical line at F — 0. As i increases, all of the 
localized (real) states move to the left until they enter either 
the left parabola of delocalized states (for genotype-1 domi- 
nated states) or the right parabola of delocalized states (for 
genotype-2 dominated states). Thus as v increases r gs will 
soon pass through at the "extinction velocity". However, 
the delocalized states do not shift to the left as v increases. 
Thus as v increases further, F g3 will eventually pass through 
r^, indicated by the center vertical line. Once this happens, 
the largest Re(r) no longer belongs to the original genotype-1 
dominated state, but rather to a genotype-2 dominated de- 
localized state. The critical velocity at which this occurs is 
the "switching velocity." In this example, since F^ < 0, the 
switching velocity is higher than the extinction velocity, and is 
therefore biological uninteresting. However, if we had > 0, 
the spectrum would look identical, except it would be shifted 
to the right. We would then have > 0, and hence there 
would always be states above Re(F) = 0. Thus there would 
be no extinction velocity, and the switching velocity would be 
biologically relevant. 



A. Solution for the Eigenvalues and Eigenfunctions 



In order to carry out the analysis sketched above, we 
must solve for the v = eigenvalues and eigenfunctions, 
and determine the values of T I and Tj . It is possible to do 
this exactly, and an outline of this calculation is presented 
in Appendix A. However, it is more straightforward and 
instructive to first examine the solutions for /ii = fi2 = 0, 
and then use perturbation theory to find the results for 
small /ii and /i2- 
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1. The /ii — fj,2 — Solution 

We first examine the system for fJ>x = fJ-2 = 0. In this 
case, the two types of individuals are completely inde- 
pendent. The problem reduces to that studied in Q. We 
use the imaginary gauge transformation to eliminate the 
velocity, and focus on the v — eigenvalues and eigen- 
states. We then have to solve the eigenvalue equation 

£ii,v— j^i.v—O ii,v—0 /i a\ 

v=o4>n = r n 4>n ■ K 1 ®) 

This problem is formally equivalent to the square well 
problem in quantum mechanics pcj . The v = right 
eigenstates are of the form 



<l>n 



where 



Gle- K 



ib 1 








(17) 



for x < -W/2 

for - W/2 < x < W/2 

for x > W/2. 



Note that the index i indicates the genotype that domi- 
nates the state. Substituting this ansatz into the eigen- 
value equation leads to 



r 



d(4, 



(18) 



We proceed by requiring that <p and its first derivative 
be continuous at x = ±W/2, which determines the con- 
stants A, B, C, and G up to an overall normalization and 
yields a transcendental equation to determine T. This 
analysis is carried out in detail in £| . The essential result 
is that provided the oasis is wide enough that a typical 
individual gives birth many times while diffusing across 
it, the ground state eigenvalue can be approximated as 
Tg^°'^ =0 ax, with a corresponding eigenfunction that 
is entirely genotype 1 and localized largely within the oa- 
sis. We will use this approximation throughout the rest 
of this paper We can also calculate the position 

of the derealization transitions T™ . These are defined 
as the amount by which the v = eigenvalue V]{ v=0 is 
shifted by the derealization velocity v™ . We thus have 

It = rjj«=° - = P>=° - £(<) 2 . Using Eq. ©, 

we find T" = hi. Note that, as claimed above, this T" is 
independent of n and depends only on i. Thus there are 
two and only two derealization threshholds, one corre- 
sponding to each genotype. 



2. Perturbation Theory in /ii and /X2 

We can now examine the results for nonzero muta- 
tion rates /ii,/i 2 > by using a non-Hermitian version 
of time-independent perturbation theory from quantum 



We 



require that [i\ and /i 2 be small, 
- < 1. The details of the 



mechanics ,20j. 
specifically that 
calculation are described in Appendix B. The eigenstates 
now all involve both genotypes, although those that be- 
gan as genotype 1 states remain dominated by this type, 
and vice versa. More precisely, we have 



4> 



l,v=0,R 



ai -A2 ~ n 



for a genotype- 1 dominated localized state, and 



,2.v=0,R 
<Pn 



(19) 



(20) 



for a genotype-2 dominated delocalized state, where ip^ 
and i/'n are gi ven m section III A II Note that we focus 
on genotype- 1 dominated localized states and genotype-2 
dominated delocalized states because no other state can 
dominate the dynamics in any regime. The eigenvalues 
are also shifted. The eigenvalues for genotype-1 domi- 
nated localized states become 



A*iai 



ai — a 2 



(21) 



while the genotype-2 dominated delocalized states have 
eigenvalues 



p2,u=0 _ p2,D=0,^t=0 



Mi a iA*2a2 
b 2 -bi 



(22) 



plus higher order terms in \x\ and /i 2 . 



B. Critical Velocities 

We can now calculate the critical velocities at which 
the dynamics changes qualitatively. There are two rele- 
vant cases. For & 2 < neither derealization transition 
occurs at positive T because neither genotype can sur- 
vive in the desert. Thus the switching velocity, though 
it exists formally, is biologically irrelevant. The extinc- 
tion velocity v e is the velocity where the ground state 
eigenvalue passes through zero. Below v e a genotype-1 
dominated population can multiply but above this ve- 
locity the population must go extinct. This velocity is 

defined by T v g % = T V = Q - & = 0, which using Eq. J2U> 
gives 



= 2y/(D ai ) 



1 - ^ 



(23) 



Note that for fix = 0, we 
the Fisher velocity for genotype 



valid to first order in \x\ and /x 2 . 
have v e — vf = 2\fUa\, 
1. 

For 6 2 > there is no extinction velocity, as genotype 2 
can survive at any velocity. However, there is a switching 
velocity v s where the population shifts from being mostly 
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Simulation Test of Analytical Result 
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FIG. 4: A comparison of the critical velocities v e and v s from 
the simulations to the analytical result, as a function of fj,i = 
10/12- The analytical result is shown as a dashed line for v s 
and a dotted line for v e . Here we have ai = 1, 52 = 0.6, 6i = 
— 1, i>2 — —0.6, and D — 0.25, where the overbars indicate the 
dimensionless units described in Appendix C. This result is 
typical of such comparisons for other parameter values. For 
these parameters, we expect the perturbation calculation of 
the analytical result to be valid for hi <C 0.4. The slight 
overestimates of v e and v s by the analytic theory for small fii 
are due to finite size effects. The underestimate as \x\ grows 
beyond the range of our perturbation expansion is due to the 
importance of higher-order terms. 



genotype 1 to mostly genotype 2. This occurs when the 
growth rate of the genotype- 1 dominated ground state 

passes through 1^ — b 2 — + . By a similar 

calculation we find 



2^D( ai - b 2 



1 



1 /LtiOi 1 jJL 2 a 2 



2 a,\ — b 2 2 a\ — b 2 



(24) 



also valid to first order in /ii and \x 2 . As mentioned above, 
this switching velocity is a sort of quasispecies transition. 
As the velocity increases passes this critical threshold, the 
population can no longer maintain the "ideal" sequence, 
even if it is below the mutational error threshold. This 
switching velocity is also well-defined for b 2 < 0, but will 
be difficult to observe because v s > v e . We can easily 
calculate the second-order corrections to both v e and v s 
(see Appendix B). 



III. SIMULATIONS 



to intuition and one check of the validity of our approx- 
imations. Second, we simulate the underlying discrete 
process, namely individuals multiplying and mutating at 
appropriate rates. This tests not just the results of our 
analysis, but also the overall applicability of continuous 
time differential equations to the real discrete popula- 
tions we model. 



A. Lattice Approximation to the Liouville 
Operator 

We made several approximations in arriving at our an- 
alytical results, including an approximations for r gs and 
for the shifts in eigenvalues with \i\ and \x 2 . To test 
these approximations, we calculate the eigenvalue spec- 
trum numerically. We discretize space, find the resulting 
discrctizcd Liouville operator, and numerically diagonal- 
ize it for a particular set of parameters. The details of 
this method are described in Appendix C. 

This approach allows us to determine the shifts in the 
growth rates as \x or v is varied, and hence the critical 
velocities v e and v s . We can compare these results to 
our analytical predictions, and therefore confirm our cal- 
culation of the critical velocities. This comparison for 
one particular set of parameters is shown in figure 0] 
Comparisons for other parameter values have been car- 
ried out, and give similar results. The eigenvalue spectra 
shown in figure |21 and figure [3] were also obtained in this 
way. 



B. Simulations of the Discrete System 

We can also simulate the underlying discrete process 
which inspired our formulation of the differential equa- 
tions we analyze in this paper. We discretize space, plac- 
ing individuals on a one-dimensional lattice which repre- 
sents our environment. These individuals move around, 
proliferate, die, and mutate. We also impose a saturation 
term so that at long times the population distribution 
settles down to a steady state. 

By comparing the steady state population profiles for 
different values of v, we determine v e or v s . These results 
can then be compared to the analytical results. This 
comparison, for one particular set of parameters, is shown 
in figure^ Note that this comparison is only for v e , as for 
these parameter values v s is not biologically relevant and 
is thus impossible to observe with this type of simulation. 
The details of this method are described in Appendix D. 



We use two different computational methods to test 
our analytic results. First, we use a lattice discretiza- 
tion of the Liouville operator C to calculate the eigen- 
value and eigenfunction spectrum for particular sets of 
parameter values. This numerical work provides an aid 



IV. DISCUSSION 

To interpret our results, it is helpful to first consider 
the behavior of a non-spatial quasispecies model. We 
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imagine a population living in a uniform environment 
whose conditions match those of the oasis. Genotype 1 
grows more quickly, but there is a mutational pressure 
away from this "ideal" genotype. Mutation away from 
this sequence is typically expected to be much more fre- 
quent than mutation back, so it is common in such mod- 
els to set [12 = 0. It is then straightforward to calcu- 
late the composition of the population. We find that the 
equilibrium ratio of genotype 1 to genotype 2 individuals 
q = Cx/c 2 is given by 



(1 - ni)ai - a 2 
A*iai 



(25) 



Thus the ratio of species 1 to 2 decreases with increasing 
\i\ until \i\ reaches the critical "error threshold" for this 
quasispecies model. At this threshold, [i\ = ai ~f 1 ' 2 , the 
"ideal" genotype can no longer survive and the popula- 
tion becomes completely dominated by genotype 2. This 
is the most famous result of Eigen's quasispecies model 

nulla. 

Our analysis finds that, in analogy to the quasispecies 
error threshold, in the spatial model there is a veloc- 
ity threshold above which genotype 1 is outcompeted by 
genotype 2. This velocity, which we call the switching 
velocity v s , is given by Eq. (|24|l . If the growth rate of 
genotype 2 in the desert is positive, we can expect to 
see this behavior in a real system. Our model naturally 
also has the traditional error threshold; for /ii > ai ~° 2 

' r ai 

and fj,2 — 0, genotype 1 is outcompeted by genotype 2 
regardless of v. It would be interesting to explore the 
interactions between the mutation-driven and velocity- 
driven quasispecies transitions. However, our analysis 
is based on the assumption of small mutation rates and 
thus focuses on the velocity-driven transition under the 
assumption that we are well away from the mutation- 
driven transition |22j| . A qualitative phase diagram of 
the different velocity-driven transitions is given in figure 

E 

Our model describes a number of important biological 
situations. The oasis could be a particularly favorable 
patch of the environment or a zone of favorable climactic 
conditions which is moving at a typical speed v due to 
seasonal weather patterns or shifts in climate. A pop- 
ulation can take advantage of this oasis by adapting to 
these conditions, but then will fare worse in the rest of 
the space. Our analysis explores how fast the oasis can 
move before the population can no longer maintain an 
adaptation to the favorable patch. Beyond this speed in- 
dividuals may occasionally find themselves in the oasis 
but cannot adapt quickly enough to benefit. 

Besides showing the existence of the velocity-driven 
quasispecies transition, our results make 5 interesting bi- 
ological predictions. We discuss each of these in turn. 

(1) The overall population growth rate decreases as v 2 
up to the critical velocity. The imaginary gauge transfor- 
mation tells us that for b 2 < 0, the exponential growth 
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FIG. 5: A phase diagram showing the transitions as a function 
of v and 62- This figure accounts for the qualitative effects of 
a nonlinear saturation term. 



rate of the population decreases as until it goes ex- 

For 62 > the exponen- 



tinct at v e = 2^a x D (l - ^ 

2 

tial growth rate of the population decreases again by 
until it reaches v s = 2^(a ± - b 2 )D (l 

Increasing the velocity further does not change the overall 
growth rate because a delocalized genotype 2 population 
then dominates. 

(2) Below the critical velocity, the genotype- 1 domi- 
nated population will extend somewhat into the desert. 
The type-1 population will extend a typical distance 
£ = — — — into the desert. From the relationship be- 
tween n and T we find that this typical distance is 

2D 



2 y /D{a l -b 1 ) 



(26) 



independent of the mutation rates. Note that this di- 



a.s 



verges 

v c ee 2^/D(ai 



(v c — v) as d 



from below, where 



bi) is the velocity at which the type-1 
dominated ground state eigenvalue r gs reaches its der- 
ealization threshold T^. This divergence will be difficult 
to observe in real biological systems because v s < v c . 

(3) The ratio of the two genotypes is proportional to 
fi. Below the switching velocity, the ratio of the number 
of type 2 individuals to type 1 approaches at long 
times. Above the switching velocity, the ratio of type 1 
to type 2 is at long times. This result, however, is 
only true within the linear model. If we add a saturation 
term to the dynamics, this term will affect the long-time 
population ratios. 

(4) The ratio of genotype 1 to genotype 2 is indepen- 
dent of v except when crossing v s . We might naively 
expect that as we increase v, the population gets driven 
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more towards the desert, and thus the population shifts 
away from genotype 1 towards genotype 2. The gauge 
transformation ensures that this does not happen. Un- 
til we cross v s , the ratio of genotypes in the ground 
state eigenfunction, and thus the population, remains 
constant. At v s there is a sharp transition and the ratio 
of the genotypes shifts radically in favor of genotype 2. 
As we continue to increase the velocity, the ratio again 
remains constant. This result, however, is also only true 
within the linear model. A nonlinear saturation term 
"softens" the transition at v s , as described in section 

GO 

(5) The "right" ratio of genotypes dominates exponen- 
tially. If the initial state of the population below the 
switching velocity is all type 2, then type 1 will take 
over exponentially with a rate equal to the difference be- 
tween T gs and the dominant type 2 eigenvalue. This 
rate is just a\ — a 2 — A'i a i + A*2<32, and is independent of 
v. Similarly, if we begin with a type-1 dominated pop- 
ulation above the switching velocity, the population will 
become dominated by genotype 2 exponentially with rate 
bi — b\ — ^,20,2 + Mi a i> again independent of v. 



A. Effects of a Nonlinearity 

Thus far we have concentrated on purely linear sys- 
tems, without much discussion of nonlinear terms which 
cause the population to saturate. This approximation 
is justified because the presence of a nonlinear satura- 
tion term will not affect the critical velocities. However, 
since all populations can be expected to saturate at some 
point, it is important to consider the general implications 
of such a nonlinearity. 

In the discussion to this point, we have assumed that 
the eigenfunction with the largest growth rate will domi- 
nate the population, so that at long times we can neglect 
all but this dominant state. In the nonlinear case this is 
still roughly true. We can account for the effects of the 
nonlinearity by using mode couplings as in Q , and antic- 
ipate that the fastest growing eigenfunction will suppress 
the others and dominate the population. There is, how- 
ever, one important exception to this idea. For the prob- 
lem considered here, the fastest-growing localized eigen- 
function vanishes exponentially outside of the oasis, and 
hence will not suppress the growth of the fastest-growing 
delocalized eigenfunction. Thus below the switching ve- 
locity v s , the population is not in fact completely dom- 
inated by genotype 1. Rather, there is a genotype-1 
dominated population inside the oasis and a genotype- 
2 dominated delocalized population in the desert. As 
we increase v we decrease the growth rate of the local- 
ized type 1 dominated state without changing the growth 
rate of the delocalized genotype 2 dominated state. Since 
the overall population sizes are set by the growth rates 
and the nonlinear terms, this will lead to a shift in the 



overall population density from type 1 to type 2. Thus 
there will be some ^-dependence in the genotype ratio 
even below v s . Just above v s , the fastest-growing de- 
localized eigenfunction will not completely suppress the 
fastest-growing localized eigenfunction, so similarly there 
will be some ^-dependence above v s . The transition at 
v s is thus softened by the presence of the nonlinearity. 
While the dominance will still shift at this critical veloc- 
ity, the transition will have some width dependent on the 
saturation term. 

The nonlinearity will also impose a maximum carrying 
capacity on the system. As the system approaches this 
carrying capacity, the growth rates will slow down. Thus 
if we start with an initial condition involving a population 
already near its carrying capacity, our analysis of the rate 
at wheh the "right" genotype composition is established 
will be an overestimate. While the results regarding the 
switching and extinction velocities and genotype compo- 
sition still hold, the dynamics of reaching the resulting 
steady states will be slower. Our results for the rates are 
based on the linearization of the nonlinear model around 
ci = C2 = 0, and so are valid as long as the ratio of the 
population to the carrying capacity is small compared to 
1. 



B. Other Geometries 

Many alternative geometries are clearly possible. One 
obvious choice is a circular oasis in an infinite desert, 
as considered in £|. The primary effect of a nontriv- 
ial two-dimensional geometry is to change the eigenvalue 
spectrum. In particular, T gs and, if the desert is not infi- 
nite, will shift. The resulting change in the switching 
and extinction velocities will depend on the growth rates 
Oj and bi and the linear size of the oasis W , but not on 
\x. However, the qualitative behavior is unaffected. Fur- 
thermore, if W is large compared to y/D/a, the typical 
length an individual diffuses before giving birth (the pre- 
cise requirement will vary with the specific geometry), 
this shift in critical velocities will vanish and the results 
will reduce to those calculated above. 



Appendix A: Exact Solution for the v — System 

It is straightforward, though tedious, to find the exact 
solution for the v = eigenfunctions and eigenvalues. We 
start with the ansatz 



4> l 



Ae~ 
Ce 



Be- 
Ee 



for ^- < x, 



<t>" = 



p e ik 1 x _j_ Q £ -ik 1 x _|_ j_f e ik 2 x _j_ j £ -ik 2 x 

Je lklX + Ke- tklX + Le ik2X + Me~ lk2X 



(27) 



(28) 
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for 



w 

2 



< X < 



, and 



,r _ i Ne KlX + Oe K2X 

9 - I p e KlX + Q e H. 2 X 



(29) 



for x < — where the 16 parameters 

A,B,C,E,F,G,H,I,J,K,L,M,N,0,P, and Q are 
constant coefficients. Demanding that this ansatz 
satisfy the eigenvalue equation yields a system of 12 
equations relating these constant coefficients, K, k, 
and the eigenvalue T. We use 8 of these equations to 
eliminate 8 of the 16 coefficients, and the remaining 4 to 
determine K\,K2,k\, and fc 2 in terms of the eigenvalue 
T. We find that for small Hi and fi 2 (tf^t* < 1 and 

(«i 



— < 1 



Dni = 



Dk\ 



Dki = - 



r — & 2 + M2«2 + 

r — b\ + niai — 

-V + ax — fixai 

r + a 2 - [l 2 a 2 



Ml a 1^2Q2 

b 2 - h 

Ml a l/ i 2Q2 

6 2 - h 

A t l a l/ i 2Q2 

ai — a 2 
_ Mi Q iM2a2 
ai — a 2 



(30) 

(31) 
(32) 

(33) 



which serves to confirm our perturbation theory calcula- 
tion. The expressions when \i\ and fx 2 are not small are 
straightforward to calculate but unwieldy to write down. 

We now demand that <fi and ^ be continuous at 



x = ± 



yielding 8 equations for the remaining 9 un- 
knowns (8 constant coefficients and T). These results 
lead to a transcendental equation for T, the solutions to 
which are the eigenvalues of the system. Choosing a par- 
ticular eigenvalue from among these possible solutions, 
we can easily determine the remaining unknowns (up to 
an overall normalization) from the rest of the equations. 
By requiring 



4> L (x) ■ cj) R (x)dx = 1, 



(34) 



we determine the normalization, and thus the exact so- 
lution. 

In practice, the transcendental equation for T is quite 
complicated and we can only solve it numerically. How- 
ever, this exact approach does provide a useful check to 
the discretized numerical solution described in section 
IIII The most important result is the fi = ground state 
eigenvalue F^ v=0 . Provided that the oasis is much wider 
than the distance a typical individual diffuses before giv- 
ing giving birth, we have T^ v=0 ss ai 0. 



Appendix B: Perturbation Theory in /ii and ^12 

We can use standard perturbation theory from quan- 
tum mechanics |2£j to determine the results for fn f [12 > 



from the fi% = fi2 = solution. We first rewrite the Li- 
ouville operator as C v= q = Cq + C±, where C\ is propor- 
tional to the (small) mutation rates \x\ and /Z2, and 



C 



_ ( Ddl + r^x) 







Dd 2 x +r 2 {x) J ' 

-/xi<zi [i 2 a 2 



'M20 2 



(35) 
(36) 



where r^x) = {ai+n i ai)9{W-\x\) + {b i -\-n i a i )6{\x\-W). 
We know the eigenvalues and eigenfunctions of Cq from 
section ITT1 and from Q. For our purposes, all that is 
necessary is that the eigenfunctions are an orthonormal 
set of left and right functions 4>„ and cj> R . 

From non-degenerate time-independent perturbation 
theory we find that the \i > eigenvalues are related 
to the /i = eigenvalues by the formula 



r„ 



+ 



(37) 



(J <t> L m {x)C^{x)dx) (J ^{x)C^{x)dx) 



plus terms of third and higher order in C\ |20j. Note 
that this result differs slightly from standard theory be- 
cause C\ is non-Hcrmitian. The first-order term in C is 
straightforward to calculate. It is simply — /^a^ for states 
dominated by genotype i. The second-order term is more 
complicated, because the \i — eigenfunctions are of the 
form 



<?r, 



1b 1 





or d) 







(38) 



where the {V^} and {V'ml each form an orthonormal set 
of eigenfunctions of the one-genotype problem. The ipn 
and ip^ are almost, but not quite, orthonormal to each 
other. 

In the approximation that these two sets of eigenfunc- 
tions are indeed orthonormal, the second order term for 
r n is easy to calculate. For genotype-1 dominated lo- 
calized states, it is MiaiA12a2 . For genotype-2 dominated 
delocalized states, it is Pl fc ° 1 _ M b 2 1 a2 ■ In the same approxi- 
mation, the corrections to the eigenfunctions are given 
by 



4> 



l.D=0,il 



^1"! ,1.2 



0,1— Q2 ^ n 

for a genotype-1 dominated localized state, and 



fc 2 -h, V n 



(39) 



(40) 



for a genotype-2 dominated delocalized state. 

These results, and the observation that T^ v=0 = a\ 
and p*^-"-° — 5 2j a llow us to calculate the critical ve- 
locities. To second order in [X\ and fi 2 , the extinction 
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velocity is given by 



v e = 2y/{D ai ) 



1 



Ml Mi 



M1M202 



(41) 



2 8 2(ai-a 2 ). 
and the switching velocity is 

Miai - M202 (Mi fl i - M202) 2 



2y/D( ai - b 2 



1 



2(oi 



8(oi 



Mi°iM2a2 



2(ai - a 2 )(ai - 6 2 ) 2(b 2 - &i)(a x - & 2 ) 



(42) 



The first order part of this result is quoted in section UTI 
above. 

This perturbation expansion relies on the assumption 
that [i\ and /7 2 are small. Specifically, we require 



< 1 



«1, 



(43) 



a 2 — ai 62 — h 

for i = l,2 for all of these results to be valid. 

Appendix C: Lattice Approximation 

Here we describe the details of the lattice approxima- 
tion to the Liouville operator. We begin by discretizing 
space, replacing it by a regular lattice with lattice con- 
stant £q. We define c" to be the population of genotype 
a at the lattice point x. We then have to solve the eigen- 
value equation 



dc a x {t) 
dt 



x' ,a>' 



C(x,a,x',a')cZ(t)=rc2(t). (44) 



The discretized version of the Liouville operator is 



£ = 



D 



' 



^2J2[e-^\x)°"*(x + l \ 



x a—1 



2D 

22 



+e— \x + e ) aa (x\ -2 cosh 
+ Y / [Ui(x)\x) 11 (x\ + U 2 (x)\x) 22 (x\ 

x 

+ /-tiai|x) 12 (x| + /x 2 a 2 |x) 91 ' 



.7; 



(45) 



where we have used the notation |x) a ^(y| to mean the 
tensor product of localized states corresponding to c" 
and eg. We define Ui(x) = (ai - max)6(W - \x\) + 
(h - vnai)9(\x\ - W), and U 2 (x) = (a 2 - (j, 2 a 2 )0(W - 
+ (b 2 — fJ. 2 a 2 )8(\x\ — W). We impose a finite size on 
the system L (with periodic boundary conditions) and a 
finite width of the oasis W, with W <C L. 

In order for the lattice approximation to be valid, the 
lattice must be fine enough that variations in the eigen- 
function <j> between lattice points is small. This means 



we must require — —j(f 



<C 1. For small v this re- 



duces to kZ£q *C 1, k™lo <C 1, and for large v we need 
<C 1 0|. These conditions are satisfied for all of the 
calculations discussed here. 

It is now straightforward to numerically solve for the 
eigenvalues and eigenstates of the lattice version of the 
Liouville operator C The results quoted here all use a 
system size L = 512^o, with an oasis width W = 10^o 
and periodic boundary conditions. 

In comparing results, it is useful to shift to dimension- 
less units. We define a dimensionless velocity 



(46) 



Note that for /ii = /j,2 = 0, the extinction velocity v e = 1. 
We then scale all the growth rates to a\ by defining 

°i 1 r b\ a 2 - b 2 - T 

01 = — = 1, 61 = — , a 2 = — , b 2 = — , I = — . 

This implies that the dimensionless diffusion coefficient 
is 



D = 



1 



(47) 



The mutation rates are already dimensionless. We use 
these redefined units in figure [21 figure El and figure 0] 



Appendix D: Discrete Simulation 

Here we describe the details of the discrete, individual- 
based simulations. We begin with a discretized spatial 
lattice containing a uniform distribution of individuals of 
both types. We also discretize time, dividing it into small 
intervals of size At. At each time, we select a lattice point 
at random. The individuals at this point can give birth, 
move due to diffusion or drift, or mutate with appropri- 
ate probabilities. The probabilities used are simply the 
coefficients of the off-diagonal elements of the discretized 
Liouville operator described in Appendix C, times At. 
We set At to be sufficiently small that the probability of 
two or more events per step is negligible. We impose a 
saturation effect (analagous to a nonlinear term in Eq. 
(@J) by setting a maximum number of individuals per 
spatial point. 

For the simulations described in figure we use the 
parameters ai = l,a2 = 0.6, 61 = — 1,6 2 = —0.6, and 
D = 0.25, where the overbars denote the dimensionless 
parameters defined in Appendix C. We use a lattice with 
512 points, with an oasis of width 10 points and periodic 
boundary conditions, and a maximum of 200 individuals 
per point. 

For these parameter values, v s is impossible to deter- 
mine because the populations are extinct at such high 
velocities. However, we can determine v e as a function of 
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and yU2- For each value of [i\ and \xi that we test, we 
plot the total number of individuals in the steady state 
distribution as a function of v. This exhibits a transition 
from some maximum number of individuals for small v 
to approximately individuals for large v. Because of 
the imposed saturation effect, the transition is not per- 
fectly sharp but rather has some small width. We define 
the transition to occur at the point at which the number 
of individuals has dropped to the number for v = 0. 



Making a different definition would shift the extinction 
velocities slightly. 

We have also run simulations for parameter values 
where v s is biologically relevant. From the w-dependence 
of the ratio of the number of individuals of type 1 to 
type 2, we can determine the value of v s in these cases. 
Although the presence of saturation broadens the transi- 
tion as in the case of v e , the results match our analytical 
predictions. 
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